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■ Interfacial roughening denotes the nonequilibrium process by which an initially flat interface 

reaches its equilibrium state, characterized by the presence of thermally excited capillary waves. 
Roughening of fluid interfaces has been first analyzed by Flekkoy and Rothman [Phys. Rev. Lett. 
£SJ . 75, 260 (1995)], where the dynamic scaling exponents in the weakly damped case in two dimensions 

were found to agree with the Kardar-Parisi-Zhang universality class. We extend this work by tak- 
ing into account also the strong-damping regime and perform extensive fluctuating hydrodynamics 



simulations in two dimensions using the Lattice Boltzmann method. We show that the dynamic 
scaling behavior is different in the weakly and strongly damped case. 



I. INTRODUCTION 

Capillary fluctuations on an interface between two fluid phases are waves that are excited by thermal noise in 
^ , the bulk [l|-|7|. From a macroscopic viewpoint, capillary fluctuations make the interface "rough" and increase the 
effective interface width. Roughening of interfaces is a fundamental aspect of nonequilibrium dynamics and has 
been widely studied in the literature (see, e.g., |8 |4l2l| for reviews). Most models for interface growth, such as the 
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Edwards- Wilkinson [131 ] or Kardar-Parisi-Zhang [14| equations, describe a purely local growth mechanism. In the case 
of an interface between two fluids, however, one can expect significant dynamical effects arising from the coupling of 
the order parameter to the hydrodynamic now field. Indeed, as has been shown in HH , the effective Langevin 
description of a roughening fluid interface is in general non-Markovian due to the surrounding flow. The roughening of 
fluid interfaces due to thermal fluctuations is potentially relevant for the stab ility of patterns that form, for instance, 
in reaction-diffusion systems [13) Oil or during phase-transitions under shear |l9l - [22l |. Also, coalescence of droplets or 
films [23| or the dynamics of wetting transitions [24[ are potentially affected by interfacial roughening. 

Thermal roughening of fluid interfaces was first studied in [l^, [l6[ based on the equations of fluctuating hydro- 
dynamics and by simulations of an immiscible lattice gas (see also [2o| for the same problem in the presence of 
surfactants). There, it was concluded that the roughening dynamics of a weakly damped interface is characterized by 
£SJ ' the scaling exponents of the Kardar-Parisi-Zhang [14| universality class. In the present work, interfacial roughening 
. of small-amplitude capillary waves is investigated in the strong-damping regime and the associated dynamic scaling 
exponents and scaling forms are derived. We show that the growth of the interfacial roughness in the strong-damping 
regime is qualitatively and quantitatively different from the weakly damped case. The theoretical predictions are 
compared against fluctuating hydrodynamics simulations of an isothermal liquid- vapor interface. 

The paper is organized as follows: In the next section, the Langevin approach to capillary fluctuations of [TBI . [l6[ 
is summarized and applied to the roughening dynamics in the strong-damping regime. Section [IIII contains results of 
fluctuating hydrodynamics simulations of a single-component two-phase fluid performed with the Lattice Boltzmann 
method. After demonstrating that both static and dynamic equilibrium properties of capillary waves are correctly 
reproduced by the simulations, the non-equilibrium roughening of a fluid interface is investigated and compared to 



the theoretical predictions. 



C3 . II. THEORY 



A. General description of capillary waves 



Capillary waves can be described in terms of a local height function h(r\\), where r 1 1 denotes a position in the 
interfacial plane (Fig. [1]). The projected area of the interface is given by L . In the two-dimensional situation we 
focus upon, r|| = x, while the perpendicular coordinate is denoted by y. In the classical capillary wave theory [l]-BJ, 
a capillary fluctuation is understood as a rigid shift of the "intrinsic" density profile. Thus, we can define the height 
function h as 



p(r) = p\nt{y- M r ||))> 



(1) 
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FIG. 1: (Color online) Capillary fluctuations of a planar interface. Sketch of the fluctuating density profile as a function of the 
lateral coordinate y. The dotted curve represents the quiescent, mean-field profile, while the dashed curve represents the fit of 
the mean-field profile to the instantaneous density profile (solid curve). The inset shows the simulation setup and coordinate 
axis. 



where p lni is the intrinsic and p(r) the instantaneous density profile. Since, in the present case, our treatment of a 
two-phase fluid is based on a Ginzburg-Landau model (see sec. HII A| . we can take as the intrinsic profile the mean-field 
solution 

Pint(y) = \{pl + Pv) + \{pL - Pv) tanh (-) , (2) 
2 2 \w / 

where p^ and py are the liquid and vapor densities, respectively, and w is the (bare) interface width. 

In our simulations, we obtain the interfacial height h by fitting p; nt to the instantaneous density profile (see Fig.[lJ. 
If the interfacial height would, instead, be determined by means of a simple crossing criterion [i.e., p(h) — (pl+ Pv) /2], 
one would pick up local density fluctuations that are present in the interface due to its finite width. These should, 
however, not be interpreted as capillary waves since they are not associated with a lateral displacement of the interface 
profile, as expressed through eq. ([T]). Local density fluctuations affect the small-scale pro perties of the interfacial 
structure and can be understood in terms of interfacial density correlation functions [26H3ll |. In this work, however, 
we shall not consider them further [70| and stick to the definition of eq. (Q~]) . 



B. Langevin theory 

For the description of non-equilibrium roughening, a Langevin approach to the capillary fluctuation dynamics is 
convenient. Such a formalism can be derived from the equations of fluctuating hydrodynamics 0, EEG3], and we 
shall base our treatment on the Langevin theory of [Tol . [lfjj , which is summarized below. The general form of the 
Langevin equation for the height fluctuations hu of an infinitely deep film turns out to be non-Markovian [TEI, [l6| , 

hu(t) = f ds X k(t - s)F k (s) , (3) 

J — oo 

where %k is a response function that, in Fourier space, is given by (l5l [l6j 

1 2ujp 

Xk w = — -— — — , with 7k (u) = — — T77T (4 

-iu;7kH + ak 2 1/2 ) 

and Fk is a random force that satisfies a fluctuation-dissipation relation 

(F k (t)FC(0)) = k B T lk (\t\) ■ (5) 

Physically, the random force F k arises from the accumulative effect of the random stress fluctuations in the bulk 
fluid up to a certain depth below the interface (cf. [7]). In the above equations, a denotes the surface tension 
and v the kinematic viscosity of the liquid. In this work, we adopt the Fourier-transform convention a(r, t) = 
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(27r) _d f dkdui exp(ik • r — iwi)ak(w), where k denotes the wavevector in the d — 1-dimensional interfacial plane (d is 
the spatial dimension) and ui is a frequency. 

In the limit of weak damping (yk — >• 0) one obtains, by expanding the square- root in 7k of eq. ((4]), 

7kwdH =-^^-2piVk^+2pi//c + 0(^ 3/2 ). (6) 
k 

Keeping in the expansion of 7k, wd only the first and the third term, the response function in the weak-damping limit 
acquires a harmonic oscillator form and is obtained from eq. Q as 

k 1 

Xk,wdM = t, 5 — ; — 77, ; ak 3 ■ ( 7 ) 

In the time domain, this results in the Langevin equation 

ak 3 

- <9 2 /i k + vk 2 d t h^ + — — h k = r k , (8) 
2p 

with rk being a Gaussian random noise source with variance (r-k(t)r^(0)) — kBT[kS'(t) + vk 3 5(t)]/ p. The contribution 
cx 5'(t) to the random force in the time-domain arises from the first term in eq. (E|) and is, therefore, present even 
in the absence of viscosity. This is one of the most characteristic features of the present Langevin theory, which, 
despite the principal harmonic oscillator form of the response function, leads to quite distinct noise-driven dynamics. 
Equation (jHJ predicts capillary waves with an oscillation frequency and a damping rate given by 

, r wd ^fc 2 . (9) 

The slight reduction of the resonance frequency u> c due to a finite damping has been neglected here. Note that, if the 
second term in eq. (j6j also is kept, the damping rate would scale cx fc 7 / 4 [Ty] . Interes ting ly, this type of scaling has also 
been derived in a few previous works based on a different theoretical approach [y, |32j. In our simulations, however, 
we observe a behavior in agreement with eq. which can be rationalized in the context of the present Langevin 
theory only if the second term in eq. (El) is neglected - as we have done above. We also remark that the relations in 
eq. (jll) agree well with experiments (33[ and other theoretical works [33, HE] ■ While this fact provides some sort of 
justification of eq. ([7]), further studies would be desirable in order to clarify the relevance of the additional terms in 
expression (E|). This, however, is out of the scope of the present work. 
In the strong damping limit, one finds for vk 2 — > oo: 

7k,dM = ±pvk - ^ + £^ + 0(,- 2 ) . (10) 

Keeping only the leading term on the r.h.s. of eq. (fTU)) . the response function follows as 

Xk,sdM = — — — n , (11) 



which results in the Langevin equation 



d t h k + ^h k = f k . (12) 
Apv 



Here, fk is a Gaussian random noise source with variance (rk(t)r^(0)) = kBT8{t)/Apvk. Equation (| 1 2|) implies a 
decay rate of 

r sd = ■ (is) 

A Langevin equation for the height fluctuations of the form of eq. (|12p has also been derived for d > 2 in [22|, |36[ . 

Weak and strong damping regimes are separated by a critical wavenumber k c . An approximate value of k c can be 
obtained by noting that for small but finite damping, the capillary wave resonance in eq. (J7J appears at a frequency 
of (w 2 — r^/4) 1 / 2 , which becomes purely imaginary if k > k c , where 

k c = ^. (14) 
pv A 
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Thus, the weak-damping regime applies to k < k c and the strong-damping regime to k > k c . Clearly, the Langevin 
equation (fT2|) - although with slightly different numerical prefactors - could also have been directly obtained by taking 
the strong-damping limit of eq. (jHJ ■ 

The dynamic correlation function C(k, t) = (hk(t)h-k(0)) follows directly from eq. ([3]). The static correlation 
function can be most easily obtained from a fluctuation-response relation as 

C(k)^(\h k \ 2 )=k B T X (k,u = 0) = ^. (15) 

This classical result of capillary wave theory can also be derived from purely geometric considerations of the energy 
cost associated with a small-amplitude interfacial distortion In fact, the above result is equivalent to a harmonic 

approximation to the interface Hamiltonian, thus describing independent capillary waves. This is a valid approxima- 
tion in the limit of small-amplitudes and large wavelengths. In the presence of gravity, a finite correlation length is 
introduced into the static structure factor, thereby cutting off the divergence at low k It is useful to remark 

that the effective capillary wave Hamiltonian can also be obtained from a spectral analysis of a Ginzburg-Landau 
type of free energy functional (26l. [28j. Such an approach has the advantage that, in principle, the density-correlation 
function in the interface can be derived from first-principles. 

Up to numerical prefactors of the order of unity, the above expressions for the oscillation frequency and damping 
agree with the results of most theories of capillary wave dynamics in the literature @, [H, Hf| |37h 39] . In these theories, 
the liquid-vapor interface is usually taken as infinitesimally thin. The expression for w c in eq. © has also been 
explicitly derived for an inviscid fluid coupled to a Ginzburg-Landau free energy functional [4C| . Also, gen eralizations 
of capillary wave theory to non-zero vapor density [HI, EJ as well as to compressible fluids [H, [39|, HH have been 
proposed in the literature. In the case of a compressible fluid, one finds that sound waves propagating parallel to the 
interface give rise to an additional resonance peak in the dynamic capillary structure factor at a frequency ui s ~ c s k. 
In the present study, this frequency is typically much larger than the resonance frequency of a capillary wave and will 
thus be neglected. 



C. Interfacial roughening 



The effective interfacial roughness is defined as the mean-square of the height amplitudes: 

w 2 (i)HIMr,i)l 2 ) = ^£<IM*)l 2 ) = / ^rr(IMi)| 2 ), (16) 

which in equilibrium in given by [H-[H HH , 

i o , , kuT I TkL , in 2D 

W* q = W\t ->■ co ~ S- x , 12 ' . „ 17 

4 cr ^log(L// ) , m 3D 

where Iq is the minimal length scale available in the system. In the 2D-case, the sum has been computed using the 
relation ^2nLi I/"- 2 — 7I ' 2 /6, whereas in 3D, it has been approximated by an integral. It is important to realize that, 
both in two and three dimensions, the interfacial roughness is diverging with the system size, although this divergence 
is very weak in 3D. The time-evolution of the interface height under the action of the random force can be computed 
directly from the Langevin eq. ([3]), explicitly assuming that i*k(s) = for s < 0, since we are interested in the growth 
from a quiescent state at t — 0. Note also that, due to causality, Xk(t) = for i < 0. 

Roughening in the weak-damping limit has already been discussed in [lH [ill], of which the essential results shall be 
summarized first. Neglecting the effect of viscosity, which is sufficient to determine the leading order behavior, the 
response function in the time-domain follows from eq. (O as [HI, [l6| 

Xk,wd(i) = sin(w c i)0(t) , (18) 
2puj c 

from which the equal-time height-correlation function can be obtained as 

(|/ik,wdW| 2 ) = / ds / ds'xk,wd(t- s)xk,wd(t~ s')(Fk{s)F£(s')) 

Jo Jo ( 19 ) 

1k B T . a 

= ^ («et) • 
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Due to the neglect of viscosity, the correlation function describes infinitely oscillating capillary waves. It is useful 



to note that, at early times (t -C l/w c 
eq. (|16|) . follows in 2D from eq. (|T9l) as H 



^k(i)| 2 ) grows esc t 2 (Fig. The time-dependent interfacial roughness, 



dk 



sin 2 (o; e f) _ 2/3 2k B T 



37r(2pcr 2 ) 1 /3 



f/.z: 



sin 2 a; 
^573" 



(20) 



where fc m ; n and fc max are the smallest and largest possible wavenumbers in the system and u> c ,mm,max are the corre- 
sponding capillary wave frequencies. In the last step, the substitution x = oj c t has been made to exhibit the leading 
time dependence. In the range where w~4ax * ^cTrain' the value of the integral is roughly constant and, con- 
sequently, W w d(t) grows like t 1 / 3 . At early times, W w d(t) grows linearly in time (Fig. [5b), reflecting the quadratic 
growth of the height-correlation function [eq. (|19p]. Note that the extent of the early-time regime decreases relative to 
the late-time regime when the size L of the interface is increased. In 3D, we have logarithmic growth, W 2 d (t) oc logt 
[l6| . The time after which the roughness reaches its equilibrium value can be estimated as a quarter of a period of 
the capillary wave with the largest wavelength, 



a/L 3 p/(16tto-) . 

In the strong- damping limit, the Fourier transform of eq. (jlip yields 

X k ,sd(*) = ^exp(-^ t )^). 

From the expression for the correlations of the random force in this limit, 

(F k (t)F£(0)} = Apvkk B T8{\t\), 
the equal-time height-correlation function results as 



(IVsdWI 2 



k B T 
ak 2 



1 — exp 



ak 



(21) 



(22) 



(23) 



(24) 



A factor of 2 has been included in the prefactor of (|24[) to recover the correct expression for the static spectrum, 
eq. (|15[) . in the long-time limit. According to eq. (|24p . the height variance grows oc t at small times until, after roughly 
a timescale of the order of the capillary time (~ rj/ak), it reaches equilibrium (Fig. [5^). In 2D, the roughness follows 
from eq. (IM1) as 



w&(t)=t 



k B T 



dx 



1 — exp(- 



= t 



k B T 



[iX-l.xJ-aT 1 ] 



x— ak n 
x—ak n 



(25) 



where a = <t/2t] and T(n,x) is the incomplete Gamma function. As long as (afc max ) _1 < t < (cifcmin) - ; the leading 
time-dependence is unaffected by the expression in the square brackets and the roughness thus grows as W s d{t) oc t 1 / 2 
until equilibrium is reached. At early times, W s d(t) is expected to grow linearly, which is related to higher-order 
frequency terms that are neglected in eq. (|11[) (see appendix for further discussion). In 3D, we find 



k B T 

4:77(7 



[-Ei(-x) + logx] 



x—ak m i n t 



(26) 



where Ei denotes the exponential integral function. In the range (afc max ) _1 <t< (afc m i n ) _1 , W 2 d (t) oc logi, which is 
preceded by a linear growth at earlier times. The roughening time in the strong-damping regime can be approximated 
as the inverse of the overdamped relaxation rate of eq. (|13[) evaluated for the largest wavenumber, 



t, 



,sd 



rjL/na . 



(27) 



Note that, at late times, the normalization of the above W 2 d (t) is slightly different from the exact equilibrium result 
based on the sum over the discrete wavemodes, eq. (flTl) . 

The origin of the different scaling behavior for the roughness in eqs. (|20[) and (f25|) lies in the different growth speeds 
of the individual height amplitudes (|/ik(i)| 2 ) m dependence of the wavenumber (Fig. H£i). In the strong-damping 
case, the growth rate scales oc 1/fc at early times, i.e., the mode with the smallest wavenumber grows fastest and the 
roughness essentially reflects the growth of this mode. In contrast, in the weak-damping case, modes of larger k grow 
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FIG. 2: (Color online) Theoretical predictions of the growth of (a) the height amplitude {\h^(t)\ 2 ) [eqs. (|19p and (|24|l ] and (b) 
the interfacial roughness W(t) [eqs. (|20|) and (|25f) ] in 2D. The solid and dashed curves correspond to the weak damping ('w.d.') 
and strong damping ('s.d.') regime, respectively. In (a), the thin curves additionally show the growth of a large wavenumber 
mode (k' > k), whose time and length is scaled with the same factors as the small wavenumber mode (thick curves) for better 
comparison. Time is scaled by the resonance frequency ui c [eq. (O] or relaxation rate F s d [eq. (|13[) ]. depending on the regime. 
In (b), time is scaled by the corresponding roughening times, eqs. (|21|l and (|27p . 



faster and the early linear growth of W is essentially due to the largest fc-mode available in the system. However, 
since large-A modes also reach equilibrium earlier [from eq. (|19p . their roughening time is oc fc -3 / 2 ], there appears a 
i^-regime in the weak-damping case where the growth of W is due to a "saturation" effect of subsequently smaller 
fc-modes, until, finally, also the mode with the smallest k has reached equilibrium. 

Neglecting the very early growth, it is seen that, in 2D, the roughness W(t) obeys a scaling relation 

W{t) = L-v,W),*Kh a = 1/2, [ Z = ^ (wCak damping) (28) 

I z = 1 (strong damping). 

The scaling index a is an equilibrium property of the roughness [eq. (|17l) ] and is independent of the damping regime. 
The scaling function w(x) reaches a constant for x — ¥ oo and behaves as 

w(x) ~ x a ' z (29) 

for small x. As noted in PHHU, in the weak-damping limit, the scaling properties of W fall into the Kardar-Parisi- 
Zhang universality class of surface growth [T^ | .l7lll However, as shown here, the Langevin theory predicts a different 
dynamical scaling index z for strong damping. [72] In three dimensions, the roughening proceeds logarithmically both 
in the weak and strong damping regime and can not be characterized by scaling exponents. 



III. SIMULATIONS 
A. Model and setup 

Our simulations are based on the fluctuating Lattice Boltzmann (LB) method [H|-[45j, which was extended to non- 
ideal fluids in [46| (see also [13, EH for further discussions). The underlying deterministic LB model, on which our 
fluctuating LB approach is based, was introduced in [49!, 50] . For general information on the LB method we refer to the 
numerous reviews available |5lT[53| . For the present purposes, it suffices to note that this method solves the equations 
of fluctuating hydrodynamics of an isothermal non-ideal fluid (40l l5^ - [58i j . The fluctuating hydrodynamic equations 
describe the evolution of a fluctuating density and velocity field p(r), u(r) and consist of a continuity equation, 

dtp = -V • (pu) (30) 

and a momentum conservation equation, 



d t (pu) + V • (puu) = 



-V-P + V-II + V-R. (31) 
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Here, II is the viscous stress tensor 

n a p = r) + dpu c 

R the random stress tensor with Gaussian correlations characterized by 
(R aP {r,t)R lS (r',t')) = 2k B T 



S(r - r')6(t - t') , 



(32) 



(33) 



and r\ and Q the shear and bulk viscosity. As a peculiarity of the LB method, 77 and £ are in general proportional to 
the density, but can otherwise be freely tuned. For generality, we have kept the spatial dimensionality d here. The 
physics of a non- ideal fluid enters the model via the pressure tensor [3I liol. l58j 



P a fs = [po - npV 2 p - -\Vp\ 2 j S a p + n{d a p)(di3p) , 



(34) 



where po is a given equation of state (bulk pressure) and k a square-gradient parameter. In fact, the above form of 
the pressure tensor P can be derived from a Ginzburg-Landau free energy functional 



dr 



Hvp| 2 + / (p) 



since 



V • P = pV 



8JF_ 
Jp~ 



(35) 



(36) 



Here, fo is a Landau free energy density, fulfilling 



Po = pdpfa - fo ■ (37) 
We take /o to be of a simple double- well form with minima at the equilibrium densities pv, Pl (3. l59l. |60| 

Mp)=13(p-Pv) 2 (p-pl) 2 - (38) 



Note that for a symmetric Ginzburg-Landau free energy functional [as in eq. (I35p ] the bending rigidity - and with it 
the Tolman correction, which describes the dependence of the surface tension on curvature [6l| - vanishes [62| ■ In the 
absence of flow and thermal noise, the equilibrium solution of eqs. (|30|) . (|3"T|) and (f38|) fulfils V ■ P = and describes, 
in the simplest case, a liquid and a vapor phase separated by a diffuse interface. The profile is given by eq. @, with 
a width of 



2 /2k 

PL- PV\ P 



(39) 



Thermal noise, which is imparted by the random stress tensor R throughout the fluid, leads to the excitation of 
capillary waves on the interface 

H HE HE HI- 0nl y 

a small number of previous simulation studies of fluctuating 
interfaces in continuum hydrod yna mic models exist, of which the ones most relevant for the present context are based 
on lattice gas automata flH [ill |25ll (see also [63|). Recently, also a fluctuating LB scheme for the Kardar-Parisi-Zhang 
equation has been introduced [641 ] - 

All our simulations are performed in two dimensions. We place a liquid stripe of size Lx x H in a rectangular 
box of size Lx x Ly, where Ly is typically ~ 2H ~ 128 lattice units (l.u.) and Lx varies between 128 and 512 l.u. 
The box is periodic in x-direction and covered by substrates at the boundaries in y-direction (see Fig. [T]). The ratio 
of the equilibrium roughness to the film thickness is less than O(10 -2 ). We expect our simulations to approximate 
the limit of infinitely deep films and have checked in a few cases, by increasing H and Ly, that our results are 
insensitive to the film height. Static properties of capillary fluctuations turn out to be quite insensitive to the specific 
simulation parameters in the present model. In contrast, regarding dynamics, satisfactory agreement between theory 
and simulation results was found to require quite large liquid- vapor density ratios (around pl/pv = 100) and intrinsic 
interface widths not larger than 5 l.u. 
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FIG. 3: Equal-time spectrum of interfacial height fluctuations on a planar one-dimensional interface obtained from LB simu- 
lations (filled circles), compared to the theoretical capillary structure factor [solid line, eq. (|15[) ]. k denotes the wavenumber in 
the plane of the interface. Simulation parameters: L = 1024, pl = 1.0, pv = 0.5, /3 = 0.11, K = 0.08, ksT = 10 -7 , r = 1.0, 
surface tension a ~ 2.7 x 10~ 3 , interface width approximately 5 l.u. 




FIG. 4: (Color online) Capillary waves on a planar interface in the weak-damping regime. In (a), the correlation function 
C(k,t) obtained from simulations (normalized to its equilibrium value) is shown for different wavenumbers k = 2irn/L, where 
n — 1 (•), 3 (■), 5 (♦), 7 (A), 9 (t). The time axis is scaled by the theoretical resonance frequency w c , eq. (9). In (b), the 
data for the capillary wave resonance frequency u> c (k) (■) and damping rate T(k) (•) are shown. The dashed lines represent 
the theoretical predictions of eq. (O, V2uj c and 2r w d, corrected by numerical prefactors. Simulation parameters: L = 128, 
p L = 1.0, p v = 0.01, P = 0.0024, k = 0.006, v = 0.00667, k B T = 10" 8 , surface tension a = 8.7 x 10" 4 , interface width ~ 5 l.u. 



B. Equilibrium properties 

Before we turn to interfacial roughening, a number of basic equilibrium properties of capillary waves are discussed. 
This not only serves as a validation of the simulation method, but is also important since, owing to the fluctuation- 
dissipation relation, equilibrium dynamics and non-equilibrium roughening are governed by the same kinetic coeffi- 
cients. In Fig. [3J the static (equal-time) capillary wave correlation function (/ik^-k) obtained from our simulations 
is shown (see also [47j|). The real-space height profile h(r) is extracted by fitting the mean-field profile, eq. ([2]), at 
each point r to the instantaneous density profile obtained from the simulation. The spectrum is computed - after 
neglecting the initial roughening period - by averaging over 2000 snapshots in a simulation running for 10 6 timestcps, 
which is around one order of magnitude larger than the largest possible relaxation time of a capillary fluctuation in 
the system, as inferred from eqs. (|2"Tj) or (l27l) (both of which give similar estimates). Perfect agreement between the 
simulation results and the theoretical capillary structure factor, eq. (1151) . is found for practically all wavenumbers. 
Note that, on a lattice, the k 2 in eq. ([L3| has to be replaced by (the negative of) the Fourier-transform of the proper 
one-dimensional discrete Laplacian, 2 — 2cos/c. The discrete nature of the Laplacian is the reason for the upturn of 
the structure factor at large wavenumbers in Fig. [3l The difference between the continuum and lattice Laplacian is 
significant only for large wavenumbers (k > 1). 

Fig. H^i shows simulation results on the capillary wave correlation function C(k, t) = (hk(t)h-u) in the weak-damping 
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FIG. 5: (Color online) Capillary waves on a planar interface in the strong-damping regime. In (a), the correlation function 
C(k,t) obtained from simulations (normalized to its equilibrium value) is shown for different wavenumbers k = 2-nn/L, where 
n — 1 (•), 3 (■), 5 (♦), 7 (a), 9 (T). The time axis is scaled by the theoretical relaxation rate, eq. (|13[) . In (b), the data for 
the damping rate in the overdamped regime are shown. The dashed line represents the theoretical prediction of eq. I|13p . 2F s d, 
corrected by a numerical prefactor. Simulation parameters: L = 128, pL = 1.0, pv = 0.01, /3 = 0.00024, k = 0.0006, v = 0.167, 
fcsT = 10 -9 , surface tension a — 8.7 x 10 -5 , interface width ~ 5 l.u. 



regime (see caption to Fig. |4] for simulation parameters) , where an oscillatory decay, 

C(k, t) = (|/i k | 2 ) exp(-r wd |t|) cos(w c i) , (40) 

with a frequency and damping rate given by eq. ©, is predicted by the theory. Since the damping rate increases 
quadratically with k, waves with larger wavenumbers are seen to oscillate for not more than a single period until 
they become practically indistinguishable from the background noise. Fig. 0J} shows the oscillation frequency and 
damping rate extracted by fitting expression (|40|) to the simulation data. The dashed lines represent \/2w c and 2r w( j, 
where the prefactors have been included here in order to conform with the expressions of standard capillary wave 
theories. It is seen that, after the correction for numerical prefactors, both quantities compare well to the theoretical 
predictions up to a wavenumber of k ~ 0.7. Deviations at higher wavenumbers are partly attributed to the fact that 
the damping is so strong that it is hard to unambiguously extract both the frequency and damping rate from the 
data. Additionally, the hydrodynamic regime, where the transport coefficients are constants, generally breaks down 
in LB at large wavenumbers jfja . l66j . 

In the case of overdamped dynamics, capillary waves decay purely exponentially: 

C(M) = (M 2 )ex P (-r sd |i|), (41) 

which is clearly seen in the logarithmic plot in Fig. [5^. Fig. [5}d shows the relaxation rate, obtained by fitting eq. (|4Tj) 
to the data. After correcting for a numerical prefactor of 2, the theoretical prediction is well reproduced up to a 
wavenumber of k ~ 0.5. Similarly to the previous case, a possible reason for the noticeable deviations at larger 
wavenumbers might be that the viscosity becomes wavenumber-dependent for larger k. It should also be noted 
that the above results are obtained in the limit of small vapor density (pv — O.Olpt). For larger vapor densities, 
the agreement between simulation and theoretical predictions is found to become worse, even when comparing to 
theoretical expressions that take into account a finite vapor density [32J, |4l| . We also observe that the discrepancies 
grow when the interfacial width is further increased. The origin of this behavior is presently unknown. Thus, in the 
future, a closer theoretical investigation of capillary wave dynamics in diffuse interface models would be interesting, 
taking also into account effects of a finite vapor density, which seems to have not been done in a sufficiently general 
way up to now [67} • 

Of course, since the above results pertain to the linear-response regime (where the fluctuation amplitudes are small 
by definition), they could have equivalently been obtained from a study of individual capillary waves in a simulation 
without thermal noise (Onsager regression hypothesis) . Previous LB studies of capillary wave dynamics made use of 
this equivalence and obtained capillary wave dispersion relations similar to the present work (4£| [6^, [6!| . In contrast 
to the equilibrium relaxation dynamics, however, the interfacial roughening phenomenon studied in the next section 
is a genuinely fluctuation induced effect. 
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FIG. 6: (Color online) Time-evolution of the interfacial roughness W in the weak-damping regime. In (a), data obtained for 
different values of surface tension, viscosity and fluctuation temperature (see Table [I]) and a fixed system size of L — 512 
are plotted. The thick dotted curve represents the theoretical prediction, eq. (|20p . Time is scaled by the roughening time 
t r [eq. (|21[) ] and the roughness is scaled by its equilibrium value W cq . In (b), simulation parameters are fixed at pl = 1.0, 
pv = 0.01, a = 8.7 x 10~ 4 , i/ = 6.7x 10~ 3 , k B T = 10~ 8 and the system size is varied as L = 32 (•), 64 (■), 96 (♦), 128 (a), 
200 (▼). The values of the scaling exponents are a — 1/2 and z = 3/2. 
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FIG. 7: (Color online) Time-evolution of the interfacial roughness W in the strong- damping regime. In (a), data obtained for 
different values of surface tension, viscosity and fluctuation temperature (see Table [J) and a fixed system size of L = 128 are 
plotted. The thick dotted line represents the theoretical prediction, eq. (|25[) . Time is scaled by the roughening time t r [eq. (|27|) ] 
and the roughness is scaled by its equilibrium value W ocl - In (b), simulation parameters are fixed at pL = 1.0, pv = 0.01, 
a = 3.6 x 10~ 5 , v = 0.17, k B T = 10~ 9 and the system size is varied as L = 32 (•), 64 (■), 96 (♦), 128 (A), 200 (▼), 300 (o), 
400 (□). The values of the scaling exponents are a = 1/2 and 2 = 1. 



C. Interfacial roughening 

Figures [6] and [7] show simulation results for the time-evolution of the interfacial roughness in the weak- and strong- 
damping regimes. It has been made sure, by choosing simulation parameters appropriately (see Table H]), that all 
wavemodes existing on the interface exclusively fall in either one of the considered regimes. In all cases, the interfacial 
width is 5 l.u. and the density ratio of pl/pv = 100 is used in order to approximate the case of zero vapor density 
(upon which the theoretical derivation is based) as closely as possible. To ensure sufficient statistical accuracy, the 
roughness is computed by averaging over 10 — 30 independent simulations. Since it was seen in Fig. [3] that the static 
capillary correlations are correctly reproduced, it is clear that the equilibrium interfacial roughness [eq. (fl"7|)] - which 
is essentially determined by the integrated structure factor - also agrees with the theoretical predictions. Therefore, 
W eq is not discussed separately here, but instead, we turn directly to the time-evolution of the roughness. 

In Fig. [Bk, the time-dependent interfacial roughness W(t) in the weak-damping case is plotted, with each curve 
corresponding to a simulation performed for different values of surface tension, viscosity and fluctuation temperature, 
keeping the system size fixed. Data collapse is achieved by rescaling time by the roughening time t r ,wd; eq. (|2ip . and 
the roughness by its expected equilibrium value W cq , eq. p7|) . We see that the overall trend of the data is correctly 
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captured by the theory (dashed curve), with a linear growth at early times and a t 1//3 -growth at late times, until 
at around a time t rw ^ the roughness attains its equilibrium value. In the crossover region between the two growth 
regimes, however, the roughness is found to grow significantly slower than predicted by the theory. In fact, the data 
seem to be more consistent with a i 1 / 4 behavior at intermediate times. This effect is more pronounced for small 
system sizes and also found to slightly depend on the chosen simulation parameters. 

In Fig. |6Jd, the roughness is shown for different system sizes between L = 32 and 200 l.u., keeping all other system 
parameters the same. By scaling time with L z (z = 3/2) and W(t) with L a (a = 1/2), it is seen that all data points 
approximately collapse onto a single master curve with a logarithmic slope of a/z = 1/3, as expected from the scaling 
form ((28)) . However, we remark that a satisfactory scaling collapse could also be achieved with an index of z = 2 
(or any other value between 3/2 and 2), corresponding to a growth behavior W ~ t 1 / 4 in Fig. [f^i- While the reason 
for these discrepancies between simulation and theory is unclear at present, we note that these deviations cannot be 
explained by the influence of a finite viscosity, since in that case one would expect to approach a t 1/ ' 2 -power-law with 
increasing viscosity and thus, find an even larger exponent than 1/3 (cf. Fig. [8]). Also, the effect seems not to be 
related to the density ratio between liquid and vapor, as we have obtained essentially the same results for different 
values of pl/pv ranging between 0.5 and 0.001. 

Fig. [7] shows the time-evolution of the interfacial roughness in the strong-damping regime. In Fig. [7^,, the data are 
obtained from simulations of varying fluid parameters but identical system size, while in [JJd only the system size is 
varied. Good agreement with the theoretical predictions [eq. (|25p ] over approximately three orders of magnitude is 
found. In particular, the data agree well with the derived scaling function (dashed curve in [7^,) and the approximate 
power-law growth of the roughness oc t 1 / 2 is recovered. At early times, some deviations from a pure power-law 
behavior are visible, which can be attributed to the neglect of terms beyond linear order in the frequency-dependence 
of the response function [see eq. (|11[) ]. This point is further discussed in appendix |A1 In contrast to the weak-damping 
case, no rescaling of time in the plot of the theoretical W(t) is found to be necessary. The scaling collapse of the data 
in Fig. [JJd is achieved with a dynamical index z = 1, confirming that the roughening dynamics in the overdamped 
case belongs to a different universality class than in the weak-damping case. 

Finally, we investigate in Fig.|S]the roughness in the crossover region from weak to strong damping. The simulation 
data (thin solid lines) in Fig. [5] have been obtained by successively increasing the viscosity from small to large 
values, keeping all other system parameters fixed. Due to the non-Markovian nature of the interface dynamics in the 
crossover regime, it is difficult to obtain an analytic expression for W(t) in this case. We observe that the simulation 
data smoothly interpolate between the overdamped (thick dotted line) and underdamped (thick dashed line) limits. 
Similarly to Fig. [5^,, deviations between simulations and theory are noticeable at intermediate times in the limit of 
weak-damping (they appear to be more pronounced here due to the smaller system size than in Fig. [6^,). 

IV. SUMMARY 

We have investigated in this work the dynamics of capillary waves and the non-equilibrium roughening of a liquid- 
vapor interface based on an effective Langevin description and by means of fluctuating hydrodynamics simulations 
using the Lattice Boltzmann method. Although roughening is a well-known mechanism of film growth, its counterpart 
in thermally excited fluid interfaces seems to have been only rarely studied [HI, [ljl HU, [25|. As a central result, we 
showed that the non-equilibrium growth of the roughness proceeds by different dynamical scaling laws, depending 
on whether the system is either weakly or strongly damped [see eq. (|2"8j) ]. In the weak-damping case, we find 
basic agreement between our Lattice-Boltzmann simulations and previous works [i~5l [l6l . |25| , which observed scaling 
exponents characteristic for the Kardar-Parisi-Zhang universality class (a = 1/2, z = 3/2). We remark, however, that 
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TABLE I: Simulation parameters (surface tension a, kinematic shear viscosity v and fluctuation temperature fegT) used in 
Figs. [6k and[7^. In both cases, pL = 1.0 and pv = 0.01. All parameters are given in l.u.. 
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FIG. 8: (Color online) Time-evolution of the roughness in the crossover region between weak and strong damping. The dashed 
and dotted curves represent the expressions for W in the weak and strong-damping limit, eq. (|20[) and eq. (|25p . respectively. 
Simulation data are shown as thin curves for visibility, corresponding to viscosities of v = 0.0033, 0.0167, 0.04, 0.067, 0.133 l.u. 
The other simulation parameters are fixed at p L = 1.0, pv = 0.01, a = 1.4 X 10 -5 , k B T = 10 -10 , L = 128. Time is rescaled by 
the roughening time corresponding to the strong-damping limit, eq. (|27[) , 
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FIG. 9: Early-time dynamics in the overdamped regime, (a) Equal-time height-correlation function, taking into account the 
leading-order frequency dependence [eq. (IA4|) ]. for two different wavenumbers k and k' > k. Time is scaled by the damping rate 
[eq. (|13|) ] corresponding to the wavenumber k. (b) Time-evolution of the roughness based on expression (|A4[) . for two different 
system sizes L and L' = 100L. Time is scaled by the crossover time r* = l/fc^ ax i/ between early- and late-time behavior. 

an dynamic exponent of z = 2 can also not be fully excluded based on the present data. In the strong damping case, 
the roughening is governed by a dynamic scaling exponent z = 1, which is characteristic for an overdamped harmonic 
oscillator driven by white noise. The scaling exponent a, related to the size dependence of the equilibrium roughness, 
is found to be equal to 1/2, in agreement with the theory. 
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Appendix A: Early-time behavior 



Relation (jlip and thus eq. (|25|) for the roughness in the case of strong damping becomes exact in the limit of infinite 
viscosity. For large but finite viscosity, eq. (|25p describes the growth only at sufficiently late times, since, by keeping 
only the leading term in the expansion of (|10[) . the high-frequency aspects of the dynamics have been neglected. In 
the case of weak damping, we expect the early-time properties to be correctly described by eq. (|20|) , since the leading 
term in the expansion of eq. ^ is already the dominant one for large frequencies. Some insights into the early-time 
growth of a strongly damped interface can be gained by retaining the leading frequency-dependent term in 7 s d- This 
results in a harmonic-oscillator form of the response function 



k 

Xk.sdM = — 
3p 



-ui 2 — ivk^kv + erfc 3 /3p 



which becomes in real space 



Xk,sd(*) = 



3pD 



exp 



-vk'H ) 2sinh I -Dt 



(Al) 



(A2) 



with D = ^/l6^ 2 fc 4 /9 — 4<rk 3 /3p. The correlations of the random force acquire now, in addition to expression fj23[) . 
a contribution proportional to the derivative of a J- function, 



(F k (t)F£(0))=k B T 



4pvk5(t) + ^rd'(t) 
k 



(A3) 



When computing the height correlation function, this term is treated as in the weak-damping case [la ]. In this way, 
we find 



k B T 
2ak 2 



1 — exp ( — k 2 vt | ( cosh Dt sinhDt 

1 1 3 J V 3D 



(A4) 



The correlation function (|A4j) is plotted in Fig. [9^, for two different wavenumbers. The early-time growth of the height- 
correlation function can be more directly assessed by neglecting the "mass term" ak 3 /3p in the response function 
JXTJ, yielding 



<IM*)I 5 



k B T 
32k 3 v 2 p 



8k 2 vt + 3 exp 



--k 2 vt 
3 



(A5) 



Thus, for small times, the height correlation function grows oc kt 2 , until for t > (vk 2 )^ 1 , crossover to linear growth, 
characteristic of the overdamped case in the infinite- viscosity limit (Brownian dynamics), occurs. As an artifact of 
neglecting the surface tension, the system roughens for an infinite time. Note that since the frequency-expansion of 
7sd has been truncated after the linear order term [eq. (I10[) ]. a different growth behavior might result at still earlier 
times. 

Figure [9Jd shows the time-dependent roughness obtained by integrating eq. (|A4[) over all wavenumbers of a finite 
system. The early-time growth of (|ft-k(£)| 2 ) is directly reflected in the linear growth of W(t), while the i^-growth 
characteristic for an overdamped system sets in after a crossover time r*. We find that t* ~ (^max^) -1 i s determined 
only by the viscosity and the largest wavenumber of the system (here, fc max = tv/Iq), but is independent of the system 
size. Thus, for v — > oo or L — > oo, the growth of the roughness will be completely dominated by the late-time behavior, 
where the scaling expressed by eq. 1(25]) holds. 
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